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Abstract 

The banded structures observed on the surfaces of the gas giants are associated with strong zonal winds alternating in direction with 
latitude. We use three-dimensional numerical simulations of compressible convection in the anelastic approximation to explore 
the properties of zonal winds in rapidly rotating spherical shells. Since the model is restricted to the electrically insulating outer 
envelope, we therefore neglect magnetic effects. 

A systematic parametric study for various density scaleheights and Rayleigh numbers allows to explore the dependence of 
convection and zonal jets on these parameters and to derive scaling laws. 

While the density stratification aff'ects the local flow amplitude and the convective scales, global quantities and zonal jets prop- 
erties remain fairly independent of the density stratification. The zonal jets are maintained by Reynolds stresses, which rely on 
the correlation between zonal and cylindrically radial flow components. The gradual loss of this correlation with increasing su- 
percriticality hampers all our simulations and explains why the additional compressional source of vorticity hardly afl'ects zonal 
flows. 

All these common features may explain why previous Boussinesq models were already successful in reproducing the morphology 
1 of zonal jets in gas giants. 

Keywords: Atmospheres dynamics, Jupiter interior, Saturn interior 
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1. Introduction 

The banded structures at the surfaces of Jupiter and Saturn 
are associated with prograde and retrograde zonal flows. In both 
planets, a large amplitude eastward equatorial jet (150 m.s"^ 
for Jupiter and 300 m.s"^ for Saturn) is flanked by multiple al- 
ternating zonal winds of weaker amplitudes (typically tens of 
m.s"^). This alternating pattern is observed up to the polar re- 
gion (Porco et al, 2003). 

Two competing types of models have tried to address the 
question of what drives the zonal winds and how deep they 
reach into the planet (see the review of Vasavada and Show- 
man, 2005). In the "weather layer" hypothesis, zonal winds are 
confined to a very thin layer near the cloud level. Such shal- 
low models, that rely on global circulation codes, successfully 
recover the observed alternating banded structures of zonal 
flows (e.g. WilHams, 1978; Cho and Polvani, 1996). While 
earlier shallow models mostly produce rather retro- than pro- 
grade equatorial jet, this has been hampered in most recent 
approaches by adding additional forcing mechanisms, such as 
moist radiative relaxation at the equator (Scott and Polvani, 
2008) or condensation of water vapor (Lian and Showman, 
2010). In the "deep models", the zonal winds are supposed to 
extend over the whole molecular envelope (-10"^ km). The di- 
rection, the amplitude and the number of jets are reproduced, 
provided thin shells are assumed and convection is strongly 
driven (Heimpel et al., 2005). 
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The in-situ measurement by the Galileo probe further stim- 
ulated the discussion. They showed that the velocity increases 
up to 170 m.s"^ over the sampling depth between 0.4 and 22 
bars (Atkinson et al., 1997). Though this result indicates that 
the amplitude of zonal winds increases well below the cloud 
level, it hardly proves the "deep models" since the probe has 
merely scratched the outermost 150 kms of Jupiter's atmo- 
sphere (Vasavada and Showman, 2005). 

Interior models of the giants planets suggest that Hydrogen 
becomes metallic at about 1.5 Mbar, corresponding roughly to 
0.85 Jupiter's radius and 0.6 Saturn's radius (Guillot, 1999; 
Guillot et al., 2004; Nettelmann et al, 2008). This has tra- 
ditionally been used to independently model the dynamics 
of the molecular layer without regarding the conducting re- 
gion. Lorentz forces and increasing density would prevent these 
strong zonal winds to penetrate significantly into the molecular 
layer, where timescales thus tend to be slower. However, shock 
waves experiments suggest that the increase of electrical con- 
ductivity is more gradual rather than abrupt (Nellis et al., 1996; 
Nellis, 2000). Liu et al. (2008) therefore argue that zonal jets 
must be confined to a thin upper layer (0.96 Jupiter's radius 
and 0.86 Saturn's radius), since deep zonal winds would al- 
ready generate important Ohmic dissipation incompatible with 
the observed surface flux. However, Glatzmaier (2008) claimed 
that the kinetic model of Liu et al. (2008) is too simplistic as 
this does not allow the magnetic field to adjust the diff'erential 
rotation, which would severely reduce Ohmic dissipation. 

Notwithstanding this open discussion, the fact that deep 
models successfully reproduce the primary properties of the 
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zonal jets speaks in their favour. These models rely on 3-D nu- 
merical simulations of rapidly rotating shells (e.g. Christensen, 
2001, 2002; Heimpel et al, 2005). When convection is strongly 
driven, the nonlinear inertial term uVu becomes influential and 
gives rises to Reynolds stresses, a statistical correlation between 
the convective flow components that allows to feed energy into 
zonal winds (e.g. Plaut et al., 2008). 

Most of these simulations use the Boussinesq approximation, 
which assumes homogeneous reference state. This seems ac- 
ceptable for terrestrial planets, but becomes rather dubious for 
gas planets, where for example the density increases by four 
orders of magnitude (Guillot, 1999; Anufriev et al., 2005). 

The anelastic approximation we adopt here, provides a more 
realistic framework, which allows to incorporate the back- 
ground density stratification, while filtering out fast acoustic 
waves (Lantz and Fan, 1999). Jones and Kuzanyan (2009) pre- 
sented 3-D anelastic simulations that showed many interesting 
diff'erences with the Boussinesq results. The zonal flows still re- 
main large-scale and deep-seated but the local small-scale con- 
vection is severely afl'ected by compressibility and strongly de- 
pends on radius. While the main equatorial band remains a ro- 
bust common feature, it seems more difficult to get stable mul- 
tiple high latitude jets, even in the most extreme parameters. 
Kaspi et al. (2009) used a modified anelastic approach neglect- 
ing viscous heating but including a more realistic equation of 
state. In their model, the zonal flows show a more pronounced 
variation in the direction of the rotation axis than in the work 
of Jones and Kuzanyan (2009). Evonuk (2008) and Glatzmaier 
et al. (2009) pointed out that compressibility adds a new vor- 
ticity source that could potentially help to generate Reynolds 
stresses that, in fine, maintain zonal flows. 

Following theses recent studies, we focus here on the efl'ects 
of compressibility on rapidly rotating convection. The main 
aim of this paper is to determine the exact influence of the den- 
sity background on the driving of zonal flows. To this end, we 
have conducted a systematic parametric study for various den- 
sity stratifications and solutions that span weakly supercritical 
to strongly nonlinear convection. 

In section 2, we introduce the anelastic model and the nu- 
merical methods. Section 3 presents the results, starting with 
convection close to the onset and then progressing into the non- 
linear regime. In section 4, we concentrate on the zonal winds 
mechanism and related scaling laws, before concluding in sec- 
tion 5. 



2. The hydrodynamical model 

2.1. The anelastic equations 

We consider thermal convection of an ideal gas in a spher- 
ical shell of outer radius r^ and inner radius r/, rotating at a 
constant frequency VI about the z axis. Being interested in the 
dynamics of the molecular region of gas giants, we assume that 
the mass is concentrated in the inner part, resulting in a gravity 
g oc 1/r^. We employ the anelastic approximation following 
Braginsky and Roberts (1995), Oilman and Glatzmaier (1981) 
and Lantz and Fan (1999). Thermodynamical quantities, such 



as density, temperature and pressure are decomposed into the 
sum of an adiabatic reference state (quantities with overbars) 
and a perturbation (primed quantities): 



p=p-\-p' ; T = T-\-r ; p = p + p\ 



(1) 



When assuming g oc 1/r^, the poly tropic and adiabatic refer- 
ence state is given by (see Jones et al., 201 1, for further details) 
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Here T and p are the background temperature and density, m is 
the polytropic index, rf = ri/ro is the aspect ratio and Np cor- 
responds to the number of density scale heights covered over 
the layer. For example, Np = 3 corresponds to a density con- 
trast of approximately 20, while the largest value explored here 
(Np = 5) implies p//po ^ 150. 

We use a dimensionless formulation, where outer boundaries 
reference values serve to non-dimensionalise density and tem- 
perature. The shell thickness d = Vo-rt is used as a length scale, 
while the viscous diff'usion time cf jv serves as the time scale, v 
being the constant kinematic viscosity. Entropy is expressed in 
units of A^-, the small imposed entropy contrast over the layer. 
The anelastic continuity equation is 



V-(pu) = 0, 
while the dimensionless momentum equation is then 



(4) 



E\ — -\-U'Vu\-\-2e,xu = -V^-F^^-^^^^-F-V-S, (5) 
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where m, p and s are velocity, pressure and entropy, respec- 
tively. S is the traceless rate-of-strain tensor with a constant 
kinematic viscosity, given by 



Jdui duj 2 
[ dxj dxi 3 
The dimensionless entropy equation then reads 



^-p\-^^^--.^ij^^^\' 



(6) 



pt{^j^+u^ Vs\ = ^ V . {pTVs) + DipS^ (7) 

where thermal diff'usivity is assumed to be constant. As stated 
by Jones and Kuzanyan (2009), in anelastic simulations with 
large density stratification, viscous heating plays a significant 
role in the global energy balance. It involves the dissipation 
parameter Di (e.g. Anufriev et al., 2005), that is given by 
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In this formulation, turbulent viscosity is expected to dominate 
over molecular viscosity and therefore thermal diff'usion relies 



mainly on entropy diffusion rather than on temperature diffu- 
sion (Braginsky and Roberts, 1995; Lantz and Fan, 1999; Jones 
and Kuzanyan, 2009; Jones et al., 2011). That is why, the tur- 
bulent heat flux has been assumed to be proportional to the en- 
tropy gradient as in the mixing-length theory (Vitense, 1953; 
Bohm-Vitense, 1958). 

In addition to the two anelastic parameters (the polytropic 
index m and the density stratification Np) and the aspect ratio 
77, the three non-dimensional parameters that control the system 
of equations (4), (5) and (7) are the Ekman number 
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the Prandtl number 



and the Rayleigh number 
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Figure 1: Radial profile of g\dsc I dr\ for different background density stratifica- 



tions iV. = [10-M,2,3,4,5]. 
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^0 being the gravity at the outer radius, Cp the specific heat, and 
K the thermal diffusivity. This definition of the Rayleigh num- 
ber is based on the entropy jump over the layer and on values 
of the physical quantities at the outer boundary. However, for 
analysing where convection sets in first, the radial dependence 
of the diffusive entropy gradient has to be considered. This can 
be cast into a depth-dependent Rayleigh number: 
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where Sc is the solution of (see Jones et al., 2011) 
V . {pfVsc) = 0. 



(12) 
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In our simulations, thermal and viscous diffusivities are as- 
sumed to be constant, but gravity and entropy vary with radius. 
Gravity is decreasing outward (as g oc 1/r^), but the entropy 
gradient is inversely proportional to p T (Eq. 13) and therefore 
increases rapidly toward the surface for large density stratifi- 
cations. As a consequence, we can see on Fig. 1 that this lo- 
cal Rayleigh number increases outward for density stratification 
larger than 3, while it increases downward for weaker density 
stratifications. 

Several authors argued that Rayleigh numbers at mid-depth 
may offer a more meaningful reference value (e.g. Unno et al., 
1960; Gough et al., 1976; Glatzmaier and Gilman, 1981). Here 
we follow Kaspi et al. (2009), who suggest to use a mass- 
weighted average defined as 
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where the average is defined as <• 
£^(---)pr^dr/ f'^pr^dr (Kaspi et al, 2009). 



(14) 



•>P 

This new 
definition of the Rayleigh number can be related to the usual 
one (Eq. 11), by the coefficient {g\ds/dr\c}p, that has an 
analytical solution for a polytropic index of 2 
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Moreover, in the limit of small viscosity and thermal diffusiv- 
ity, the force balance is expected to become independent on vis- 
cosity (low Ekman number). Therefore, following Christensen 
(2002) and Christensen and Aubert (2006), it is more relevant 
to consider a modified Rayleigh number that does not depend 
anymore on diffusivities, defined as 
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To get a non-dimensional measure of the radial heat transport 
that is also independent of the thermal diffusivity, we use a mod- 
ified Nusselt number 



Nu* = 



pTQd^(dSc/dr) 



= NuEPr" 



(17) 



where q is the heat flux. Once the simulation has reached the 
nonlinear saturation, the time-average of this modified Nusselt 
number should become constant with radius as there is no in- 
ternal heat source in the system. Therefore, there is no need to 
consider mass-weighted average of this parameter. 

Finally, it is also useful to analyse our simulations in terms 
of an alternate Rayleigh number based on the heat flux rather 
than on the entropy contrast As 
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This flux-based Rayleigh number can be related to the modified 
Nusselt number as follows 
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dr 



Nu* = Ra* Nu* 



(19) 



Once again, this number varies with radius, and therefore it 
is useful to consider the mass-weighted average counterpart of 
this quantity 



<Ra;>p = {Ra*)p Nu*. 



(20) 



2.2. The numerical method 

The numerical simulations of this study have been carried out 
with a modified version of the code MagIC (Wicht, 2002). The 
new anelastic version has been tested and validated against dif- 
ferent compressible convection and dynamo benchmarks (Jones 
et al., 201 1). To solve the system of equations (4), (5) and (7) in 
spherical coordinates (r, Q, 0), the mass flux pu is decomposed 
into a poloidal and a toroidal contribution 



Table 1 : Values of the critical Rayleigh numbers (and flux-based counterparts) 
and azimuthal wavenumbers for diff'erent background density stratifications 
(withE= 10-4 and Pr= 1). 



A^n 



Raci 



Ra! 



<Ra;)p 



10-2 1.739 xlO^ 

1 5.175 xlO^ 

2 1.141 xlO^ 

3 1.529 xlO^ 

4 1.852 xlO^ 

5 2.341 X 10^ 



1.739x10-^ 
5.175 xlO-'7 
1.141 X 10-6 
1.529x10-6 
1.852x10-6 
2.341 X 10-6 



2.666 X 10-^ 21 

8.275 X 10-^ 34 

1.496x10-6 53 

1.326x10-6 72 

9.023 X 10-^ 83 

5.727 X 10-^ 93 



pM = V X (V X Wer) -\-V xZer 



(21) 



where W and Z are the poloidal and toroidal potentials. W, Z, p 
and s are then expanded in spherical harmonic functions up to 
degree ^max in the angular variables 6 and and in Chebyshev 
polynomials up to degree Nr in the radial direction. The equa- 
tions for W and Z are obtained by taking the horizontal part of 
the divergence and the radial component of the curl of Eq. (5), 
respectively. A more exhaustive description of the numerical 
method and spectral transforms involved in the computation can 
be found in (Gilman and Glatzmaier, 1981) and in (Christensen 
and Wicht, 2007). 

2.3. Boundary conditions 

In all the simulations presented in this study, we have as- 
sumed constant entropy and stress-free boundary conditions for 
the velocity at both limits. Under the anelastic approximation, 
the non-penetrating stress-free boundary conditions is modified 
compared to the usual Boussinesq one as it now involves the 
background density 



d I 1 dW\ _ d I 1 
dr \ r^p dr I dr\ r^p 



for r= {rural (22) 



The non-penetrating condition of the radial motion is not af- 
fected by the anelastic approximation and remains t/^ = (or 
ly = 0) at both limits. 

3. Physical properties of compressible convection 

We have fixed the aspect ratio // = 0.6 to the lower value 
suggested for the molecular to metallic hydrogen transition in 
Jupiter (0.85 Rj) and Saturn (0.6 Rs ) (e.g. Liu et al., 2008). This 
limits the numerical difficulties associated with thinner shells 
and therefore allows to reach higher density stratification. The 
Ekman number is kept fixed at E = 10-^, which is larger than 
the most extreme values chosen by Jones and Kuzanyan (2009) 
but allows us to conduct a large number of simulations. Fol- 
lowing previous Boussinesq studies (e.g. Christensen, 2002), 
the Prandtl number is set to 1, while the polytropic index is 
m = 2 following Jones and Kuzanyan (2009). Being mostly in- 
terested in the effects of density stratification on zonal flow gen- 
eration, we have performed simulations at five diff'erent values 
of Np (10-^, 1, 2, 3, 4, 5) that span the range from nearly Boussi- 
nesq to a density contrast of 150. In the gas giants, the density 



jump between the 1 bar level and the bottom of the molecu- 
lar region corresponds to Np ^ 8.5 (e.g Guillot, 1999; Nettel- 
mann et al., 2008). For numerical reasons, we cannot aff'ord 
to increase Np beyond 5. However, since the density gradient 
decreases rapidly with depth, Np = 5 covers already 99% of 
the molecular envelope when starting at the metallic transition. 
The remaining rather thin upper atmosphere (roughly 500 kms) 
may involve additional physical eff'ects such as radiative trans- 
fer, weather eff'ects or insolation, that we do not include in our 
model. For each of the stratifications used in this study, the 
Rayleigh number has been varied, starting from close to on- 
set up to 140 times the critical value. Altogether more than 
140 simulations have been computed, each running for at least 
0.2 viscous time ensuring that the nonlinear saturation has been 
reached (see Table 2). 

The numerical resolution increases from (A^^ = 73, ^max = 
85) for Boussinesq runs close to onset to (A^^ = 161, ^max = 
341) for the more demanding highly supercritical simulations 
with strong stratification. In the latter cases, we have resolved 
to imposing a two-fold or a four-fold azimuthal symmetry, ef- 
fectively solving for only half or a quarter of the spherical shell, 
respectively. A comparison of testruns with or without symme- 
tries showed no significant statistical diff'erences. 

3. 1 . Onset of convection 

The linear stability analysis by Jones et al. (2009) shows that 
density stratification can have a significant influence both on 
the value of the critical Rayleigh number and on the location of 
the first unstable mode. To facilitate the comparison between 
diff'erent density contrasts, we have therefore computed the cor- 
responding critical Rayleigh numbers (Tab. 1). The values have 
been obtained by trial and error, monitoring the growth or decay 
of initial disturbances at various Rayleigh numbers. 

Figure 2 emphasises how the location of the convective insta- 
bility changes when Np is increased. In the nearly Boussinesq 
case {Np = 10-^), convection sets in near the tangent cylinder. 
The convective columns show a strong prograde spiralling, typ- 
ical for rapidly rotating convection at a moderate Prandtl num- 
ber (e.g. Zhang, 1992). For Np = I, the wavenumber increases 
but the instability stays attached to the inner boundary. How- 
ever, when increasing Np further, it moves outward and is pro- 
gressively confined to the very thin outermost region. 

This behaviour can be explained by the depth-dependent 
of buoyancy expressed via the modified Rayleigh number 
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Figure 2: Radial component of the velocity Ur in the equatorial plane for different density stratifications very close to onset of convection (less than 5% over the 
critical Rayleigh number). Outward flows are rendered in red, inward flows in blue, dimensionless radial velocities are expressed in terms of Reynolds numbers. 



(Eq. 12). Figure 2 demonstrates how larger buoyancy values 
are more and more confined close to the outer boundary as Np 
grows. This implies a decrease in radial lengthscale, which 
goes along with a rapid increase in wavenumber. This allows 
the instability to maintain a more or less spherical cross section 
(Glatzmaier and Gilman, 1981; Jones et al., 2009). 

3.2. Convection regimes 

The properties of rotating thermal convection have been ex- 
plored extensively under the Boussinesq approximation (e.g. 
Zhang, 1992; Grote and Busse, 2001; Christensen, 2001). The 
respective studies show that when the Rayleigh number is in- 
creased beyond onset distinct regimes with diff'erent tempo- 
ral behavior are encountered. Here we explore whether these 
regimes are retained in the presence of density stratification. 

Figure 3 illustrates the time dependence in the temporal vari- 
ation characteristic for the diff'erent regimes at A^p = 2 and in 
Fig. 4 we attempt to show how the regimes boundaries change 
on increasing Np. At Np = 2, we still find the same regimes 
identified in the Boussinesq simulations. Very close to onset, 
the solutions simply drift in azimuth due to its Rossby wave 
character (Busse, 1994) and the kinetic energy becomes sta- 
tionary (Fig. 3a). The drift is a consequence of both the height 
change of the container and the background density stratifica- 
tion encountered by the convective instability. Its amplitude and 
direction depends on the chosen stratification (Evonuk, 2008; 
Glatzmaier et al., 2009). At a still rather small supercriticality, 
the solution starts to vacillate in amplitude, while retaining its 
columnar form. Figure 3b demonstrates that the poloidal and 
toroidal kinetic energies oscillate in phase and have a compara- 
ble magnitude. 




(Ra-Ra,,,)/Ra„,, 

Figure 4: Regime diagram (Ra,A/p) displaying how the difl'erent nonlinear 
regimes of rotating compressible convection depend on supercriticality and 
density stratification. Each symbol corresponds to one simulation computed 
in this study. The five black triangles denote the localisations of the five simu- 
lations emphasised in Fig. 3. 
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Figure 3: Time series of the kinetic energy for different simulations with a moderate density contrast {Np - T). In each panel, the dotted black line corresponds 
to the total kinetic energy, while solid black line and solid gray line correspond respectively to toroidal an poloidal parts of the kinetic energy. The simulation 
displayed in (a) corresponds to Ra = 1.145 x 10^ (1.004x supercritical), (b) to Ra = 1.2 x 10^ (1.05x supercritical), (c) to Ra = 2.2 x 10^ (1.93x supercritical), (d) 
to Ra = 8 X 10^ (7x supercritical) and (e) to Ra = 2 x 10^ (17.5x supercritical). These five simulations can be localized in the regime diagram displayed in Fig. 4. 



A further increase of the Rayleigh number leads to a third 
regime with a more chaotic time-dependence (Fig. 3c). The 
zonal flow production is already so efficient that the toroidal 
energy dominates the kinetic energy budget. Since the poloidal 
flow is directly driven by convective up- and downwellings, its 
amplitude can serve as a proxy for the convective flow vigor. 
Zonal flows are the axisymmetric contribution of toroidal flow 
and already dominate the toroidal energy here. Both the con- 
vective and the zonal flows display large amplitude modula- 
tions. A fourth regime with intermittent nearly oscillatory con- 
vection is encountered once Ra > 4 x Ra^ (panel (d) in Fig. 3). 
Here, the zonal flow periodically becomes so large that the as- 
sociated shear disrupts the convective structures (see for fur- 
ther details Christensen, 2001; Grote and Busse, 2001; Simitev 
and Busse, 2003; Heimpel and Aumou, 2012). The Reynolds 
stresses then decrease and the zonal flow cannot be maintained 
any more against viscous forces. The convection recovers once 
the zonal flows amplitude has decreased sufficiently. A simi- 
lar but less regular competition may also cause the large varia- 
tion in the chaotic regime. Typical for this intermittent regime, 
which has mainly been observed in Boussinesq simulations but 
has also been reported in some anelastic simulations of turbu- 
lent convection in solar-like stars (Ballot et al., 2007), is that the 
convective features typically occupy only an azimuthal fraction 
of the volume. Finally, in the strongly supercritical regime, con- 
vection becomes once more chaotic but the variations are much 
smaller than in the first chaotic regime (Fig. 3e). Convective 



features now fill the whole spherical shell and efficiently drive 
strong zonal winds. 

For A/p = 2 we thus find the following regime succession also 
typical for Boussinesq convection at Pr < 1 (Grote and Busse, 
2001): 

Drifting -^ Vacillating -^ Chaotic -^ Intermittent -^ Chaotic . 

The regime diagram displayed in Fig. 4 shows that the first 
chaotic and the intermittent regime tend to vanish for larger val- 
ues of A/p. For A/p = [4 - 5], we found neither the first chaotic 
nor the intermittent regime, which could mean that both have 
retreated to a rather small region of the parameter space or are 
missing altogether. Due to the lack of an intermittent regime, 
however, the separation of the first and second chaotic regime is 
less clear cut, being based only on the variation amplitude. The 
narrowing of the regime diagram for increasing density strat- 
ification was also found in 3-D simulations of convection in 
rapidly rotating isothermal spherical shells (Tortorella, 2005). 

The large amplitude oscillations in the first chaotic and the 
intermittent regime require a more global interaction between 
convection and zonal flow. When the density stratification is 
increased, however, the solution becomes more small scale and 
chaotic and the convective columns progressively loose their 
integrity in the direction of the rotation axis (Glatzmaier et al., 
2009; Jones and Kuzanyan, 2009). This likely prevents the re- 
quired more global interaction and the related regimes at larger 
density stratification. 



3.3. Flow properties 

The background density stratification causes significant 
changes in the convective flow that we discuss first before de- 
scribing the zonal flows in section 3.3.2. 
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Figure 5: Radial velocity Ur on spherical surfaces at three different radii close 
to the inner boundary (r = 0.65ro), at the layer middle, and close to the outer 
boundary (r = 0.95ro). Outward flows are rendered in red, inward flows in blue, 
dimensionless radial velocities are expressed in terms of Reynolds numbers. 
From a model with No = 5 and Ra = 8 x 10^. 



3.3.1. Convective flows 

Figure 5 iUustrates how the radial velocity varies with depth 
in a strongly- stratified anelastic simulation at A/p = 5. The 
Rayleigh number of Ra = 8 x 10^ guarantees that the solution 
is safely located in the second chaotic regime. Deeper within 
the interior, the convection is characterized by larger azimuthal 
length scales and lower amplitudes. With increasing radius, 
the azimuthal length scale decreases, while the amplitude in- 
creases. The convective features remain roughly aligned with 
the rotation axis at any radius. The changes in length scale are 
coupled to the background density profile as the density scale 
height defined by Hp = -(dlnp/dr)~^ is a good proxy for the 
typical size of the convective eddies. At A/p = 5, this scale 
height decreases by a factor of five when going from the inner 
to the outer boundary, which about explains the variation shown 
in Fig. 5. 

The increase in amplitude can be explained by the varia- 
tion of the background density: since p decreases by a factor 
of about 150 from r^ to r^, a rising (falling) convective plume 
compensates this with increasing (decreasing) flow amplitude. 
Figure 5 shows an increase of about 60 which seems somewhat 
on the low side. However, the radial increase in the number of 
convective columns and the scale change show that the convec- 
tive features hardly reach through the whole shell without being 
modified. 

3.3.2. Zonal flows 

Evonuk (2008) and Glatzmaier et al. (2009) argue that the 
density stratification can significantly influence the generation 
of zonal winds. Exploring the possible related efl'ects is the 
main focus of our study. Figure 6 shows how the zonal flows 
evolve on increasing the Rayleigh number in moderately strat- 
ified simulations (A/p = 3). Close to the onset of convection, 
a thin prograde equatorial jet develops in the very outer part 
of the shell. The latitudinal extent of this jet is smaller than 
in Boussinesq convection because the convective columns are 
confined closer to the outer boundary (see the simulation with 
A/p = 3 in Fig. 2). It is flanked by weak amplitude retrograde 
jets that strongly vary along z. As Ra is increased, convec- 
tion starts to fill the whole shell, the equatorial jet broadens 
and the zonal winds become more and more geostrophic with 
a second retrograde jet that reaches throughout the shell. For 
Ra = 5 X 10^ (3.3x critical) the equatorial jet has reached a 
maximum width that it retains at even higher Rayleigh num- 
bers. The retrograde jet broadens toward higher latitudes until 
additional flanking prograde jets develop below and above the 
inner boundary at large Rayleigh numbers and confine it to mid 
latitudes. The amplitude of these high-latitude jets becomes 
comparable to that of the equatorial jet for strongly nonlinear 
simulations (Ra = 4 x 10^, 26 x Ra^), while the retrograde jet 
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Figure 6: Zonally averaged azimuthal velocity in the meridian plane for simulations with Np = 3 and different Rayleigh numbers. Colorscales are centered around 
zero: prograde jets are rendered in red, retrograde jets in blue. Prograde contours have been truncated in amplitude to emphasize the structure of the retrograde 
flows. Dimensionless azimuthal velocities are expressed here in terms of Rossby numbers. The extrema of the zonal flow velocity are indicated in the center of each 
panel. 
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Figure 7: Zonally averaged azimuthal velocity in the meridian plane for simulations with difl'erent density stratifications. 
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vigor intensifies to about 50% of this ampHtude. In this multi- 
ple jets regime, retrograde zonal winds are anchored to the tan- 
gent cylinder in a similar way to what was already observed in 
Boussinesq simulations (e.g. Christensen, 2001; Heimpel et al., 
2005). 

The strong z-dependence at low Rayleigh numbers is a typ- 
ical thermal wind feature, which can be understood by consid- 
ering the 0-component of the curl of the momentum equation 
(5): 






= rsm 



\rsm6/ 



\rsm6) 
Ra rl ds 



2 dU(p 



[Vx(p-V.S)l 



(23) 



Here D/Dt corresponds to the substantial time derivative and 
oj(P is the azimuthal vorticity component. For low Rayleigh 
numbers, inertial contributions can be neglected, and the z- 
dependence of the zonal flows is ruled by the thermal wind bal- 
ance: 

2^^Ra^^ 

E dz Pr r3 de' ^ ^ 

where overlines denote an azimuthal average. This allows us to 
define the typical lengthscale 



dz 



U(p 



Ra E rl d-s 



(25) 



For vanishing Ekman numbers, L^ approaches infinity and zonal 
flows become independent of z, recovering Taylor-Proudman 
theorem and therefore a geostrophic structure. For the low- 
est Rayleigh number solution shown in Fig. 6a, the thermal 
winds caused z- variation is still clearly apparent where t/^ is 
weak and L^ therefore small. At larger Rayleigh numbers, the 
thermal wind balance still holds in a statistical sense. Accord- 
ing to Eq. (24), du^/dz increases linearly with Ra, while, at 
least closer to onset, the zonal flow amplitude rises faster. L^ 
therefore grows with the Rayleigh number and the z- variation 
is gradually lost. In the simulations of Kaspi et al. (2009) how- 
ever, the zonal winds still vary strongly close to the surface. 
These are likely promoted by the depth-dependence of the ther- 
mal expansion coefficient, which is neglected in our ideal gas 
model. 

Figure 7 compares the zonal flows at different density strat- 
ifications for simulations with similar <Ra*)p and demonstrates 
that the structure changes very little for A/p < 3. In the two cases 
at A^p =4 and A/p = 5 the high-latitude jets are still little pro- 
nounced and the retrograde jet is also rather weak. According 
to the development shown in Fig. 6 they seem to correspond 
to a lower supercriticality. The fact that A/p = 4 and 5 simu- 
lations have a supercriticality that is about three times higher 
than that of the Np = 2 case, however, implies that significantly 
larger supercriticalities are required to reach a comparable mul- 
tiple jet structure when the density stratification is strong. Since 
both large Rayleigh numbers and strong density stratifications 
promote smaller scales, we could not afford to further increase 
Ra for A/p > 3. While the zonal winds are very geostrophic for 
smaller density stratification, additional non-geostrophic small 



scale features start to appear for larger stratifications as it has 
already been reported by Jones and Kuzanyan (2009). They 
are strongly time-dependent and therefore disappear when time 
averaged quantities are considered. 

Except for these secondary features and the requirement for 
larger supercriticalities at stronger stratifications, the zonal flow 
structure is surprisingly independent of the background density 
profile. We now turn to discussing their amplitude. 

4. Scaling laws 

To evaluate how the different flow properties depend on the 
background density stratification, we consider time averaged 
quantities in the following. Each simulation has been run for 
at least 0.2 viscous diffusions times to get rid of any transients 
and to allow for an averaging period long enough to suppress 
short term variations. The Reynolds numbers Re and Rossby 
numbers Ro = E Re, used to quantify the amplitude of differ- 
ent flow contributions are based on the non-dimensional time 
averaged root-mean square velocity: 



with 
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(26) 
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where r is the averaging interval, Vs is the volume of the spher- 
ical shell and dV is a volume element. As pointed out by Kaspi 
et al. (2009), a comparison of simulations with different den- 
sity backgrounds may be more relevant if one considers mass- 
weighted average quantities. Therefore, the mass- weighted 
counterparts of Reynolds and Rossby numbers are defined us- 
ing the time averaged kinetic energy, that is 
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(28) 



(29) 



Figures 8 and 9 show how diff'erent time-average proper- 
ties for the various density stratifications explored here depend 
on the flux-based Rayleigh number (Eq. 19) and the mass- 
weighted flux-based Rayleigh number (Eq. 20), respectively. 
The total Rossby number Ro (Fig. 8a) and the poloidal Rossby 
number Ropoi (Fig. 8b), as well as their mass-weighted coun- 
terparts (Fig. 9a-b) first rise steeply. The slope flattens once 
convection fills the whole shell and seems to reach an asymp- 
totic value for large Rayleigh numbers. The asymptotic slopes 
are larger for the poloidal than for the total flow amplitude. 
In both Figs. 8 and 9, the curves collapse at larger Rayleigh 
numbers, and the asymptotic slopes become independent of 
the density stratification. This is even more convincing when 
mass-weighted average quantities are considered, which basi- 
cally confirms Kaspi et al. (2009). The remaining differences 
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Figure 8: Time-averaged properties plotted against flux-based Rayleigh number defined in Eq. (19): (a) Rossby number; (b) Rossby number for the poloidal flow 
component; (c) ratio of square velocity over non-zonal square velocity; (d) modified Nusselt number. 
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Figure 9: Time-averaged properties plotted against mass-weighted flux-based Rayleigh number defined in Eq. (20): (a) mass-weighted Rossby number; (b) mass- 
weighted Rossby number for the poloidal flow component; (c) ratio of kinetic energy over non-zonal kinetic energy; (d) modified Nusselt number. 
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at low Rayleigh numbers can at least partly be explained by the 
different critical values. Compared to the gas giants, the Rossby 
numbers where our simulations reach the asymptotic regime are 
roughly one order of magnitude larger. However this is likely a 
consequence of the moderate Ekman number used in our mod- 
els as Christensen (2002) showed that this transitional Rossby 
number decreases with Ekman number. 

Figure 9c displays the ratio of the total to the non-zonal ki- 
netic energy, while Fig. 8c corresponds to the ratio of square 
velocities only. For all density stratifications, this ratio initially 
increases in the weakly supercritical regime but decays once the 
convection becomes strongly nonlinear. Similarly to what has 
been reported for Boussinesq simulations (Christensen, 2002), 
the maximum is reached close to the transition to the intermit- 
tent regime (see Fig. 4 for A/p < 3). 

In the nearly Boussinesq simulations, the maximum ratio is 
Eydn/Enz - 20, which means that 95% of the energy are carried 
by zonal flows. For Np = 5 this decreases to about 90% at a 
maximum ratio of E^^/E^j^ ^10. These values are on the low 
side of the expected ratio for Jupiter: around 50 - 100 from 
the work of Salyk et al. (2006). However, this moderate maxi- 
mum ratio is once more a consequence of the too large Ekman 
number used in our simulations (Christensen, 2002). The slope 
of the decrease for larger Rayleigh numbers also seems to vary 
with A/p, but a clear dependence is difficult to extract from our 
dataset. Considering mass-weighted quantities slightly helps to 
merge the different simulations, but they are far from collapsing 
as nicely as for the Rossby numbers. 

Panels (d) in Figs. 8 and 9 show the dependence of the modi- 
fied Nusselt number on the Rayleigh numbers. Close to onset of 
convection, the heat transport is dominated by conduction. The 
classically defined Nusselt number remains close to unity and 
the Nusselt number increases only weakly with the Rayleigh 
number. When the convective motions become more vigor- 
ous, convective heat transport starts to dominate and the Nusselt 
number rises steeply. Using the mass-weighted Rayleigh num- 
ber <Ra*)p allows to collapse the curves for different density 
stratifications and a similar asymptotic slope is reached. This 
once more emphasises that mass-weighted quantities are prob- 
ably more meaningful parameters, when considering strongly 
nonlinear convection. 

Based on the results displayed in Fig. 9 we have estimated the 
asymptotic dependence of (Ro)p, (Ropoi)p and Nu* on (Ra*)p 
to: 

(Ro)p = 2.45((Ra;)p)'^', (30) 
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Figure 10: Time average of Reynolds stresses (solid grey line) and viscous force 
(dashed black line) integrated over cylinders (Eq. 34) for a model with Np = 3 
and Ra = 2 X 10^. The vertical line corresponds to the tangent cylinder. 



The amplitude of the convective flow, quantified via the 
poloidal Rossby number, is thus found to depend on the mass- 
weighted flux-based Rayleigh number with a slope of 1/2. This 
is compatible with both laboratory experiments of turbulent 
rotating convection (e.g. Fernando et al., 1991) and previous 
anelastic numerical simulations (Showman et al., 2011) but 
slightly difl'ers from the exponent 2/5 suggested by Christensen 
(2002). Investigating the reason for this difl'erence would prob- 
ably require more simulations at lower Ekman numbers and lies 
beyond the scope of this study. 

The Nusselt number scaling Eq. (32) is in a very good agree- 
ment with previous laws derived by Christensen (2002) and 
Christensen and Aubert (2006) for Boussinesq simulations. In 
the highly supercritical regime, the convective transport domi- 
nates difl'usion. Then Nu ~ Ucom ^s, or Nu* ~ Ropoi Ss where 
Ss is a measure for the typical entropy fluctuations. Since Nu* 
and ROpoi depend with a very similar exponent on <Ra*)p, the 
typical entropy fluctuations Ss may vary very little with <Ra*)p 
in the asymptotic regime, in agreement with previous Boussi- 
nesq results (Christensen, 2002). 

In the highly supercritical regime, the majority of the kinetic 
energy is carried by zonal flows. The asymptotic law for the 
total Rossby number (Eq. 30) therefore applies to the amplitude 
of zonal winds. Based on energy considerations for the highly 
supercritical regime. Showman et al. (2011) derived a scaling 
law for strong zonal jets that suggest an exponent of 1/2 rather 
than the 1/3 promoted here. We further discuss this difl'erence 
in the following section. 



Nu* = 0.086 ((Rapp)' 



0.53 



(32) 



Black lines in Fig. 9 show the respective asymptotes. The same 
slopes have been obtained for the usual Rossby and poloidal 
Rossby numbers in Fig. 8, as the nearly Boussinesq simula- 
tions hardly depend on the density background (therefore mass- 
weighted counterparts of Rossby and poloidal Rossby numbers 
are equivalent to their usual definitions). 



4.1. Force balance and correlation of the convective flow 

To explore how the strong zonal flows are generated and 
maintained, we analyze the azimuthal component of the Navier- 
Stokes equation (5) integrated over cylindrical surfaces of ra- 
dius s. The Coriolis term is proportional to the net mass flux 
across these surfaces and therefore vanishes. The integrated 
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Figure 11: Reynolds number of axisymmetric zonal flows plotted against 
poloidal Reynolds number. 



Figure 12: Correlation C^^ plotted against mass-averaged flux-based Rayleigh 
number. The solid line is proportional to (Ra* V 
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Navier-Stokes equation thus simplifies to 



I 



P^-sd(pdz 
s ot 



Js s^ ds 



pUsU(p sd(pdz 



2- ^ (^<P] 



(33) 



sd^dz- 



On time average the two terms on the right hand side, Reynolds 
stress and viscous stress, should balance. This leads to the fol- 
lowing condition 

\ ps-—[ — \dz= \ pUsU^dz, (34) 

where the overlines indicate an azimuthal average. Figure 10 
compares the left and right hand side of Eq. (34) for a case with 
Np = 3 and Ra = 2x10^ (corresponding to <7?a*)p = 1.7 x lO'^) 
that lies very close to the asymptotic regimes found in Fig. 9. 
These two terms balance to a good degree, which proves that 
zonal flows are indeed maintained by Reynolds stresses against 
viscous dissipation. 

According to Eq. (10), we can also directly relate the viscous 
contribution in Fig. 9 to zonal flow gradients. Zonal flows are 
thus expected to increase outside the tangent cylinder, located 
at ^ = 1.5, and to mainly decrease inside the tangent cylinder. 
A conversion to Reynolds stresses suggests that the strongest 
stresses can be found around the tangent cylinder and close to 
the outer boundary. The stresses are negative around the tangent 
cylinder and predominantly positive elsewhere. All this corre- 
sponds nicely to the multiple jet cases described above with 
a prograde equatorial jet, a retrograde jet around the tangent 
cylinder, and a prograde, high latitude jet in each hemisphere. 

Since the typical length scale of the zonal jets remains of or- 
der unity, Eq. (10) suggests the scaling (e.g. Christensen, 2002) 



(Rezon)p ~ Cs(f) (RCpolV 



(35) 



where (Rezon)p is the Reynolds number of axisymmetric zonal 
flows and C^^ quantifies the correlation of Us and U(p throughout 
the shell: 
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Here, the rectangular brackets [• • • ] corresponds to an average 
over s and z. In Fig. 11 we show (Rezon)p versus (Repoi)p for 
all the simulations of this study, while Fig. 12 shows the depen- 
dence of Cs(p on the mass-weighted flux-based Rayleigh number 
for some snapshots. 

For lower Rayleigh numbers and thus smaller poloidal 
Reynolds numbers (Rezon)p scales roughly like (Repoi)p and C^^ 
is largely independent of the Rayleigh number. Convection as- 
sumes the form of regular large scale columns that are tilted in 
prograde direction due to the height change in the container and 
the background density stratification (e.g. Busse, 1983, 2002). 
The integrity of the rolls in z-direction and the consistent tilt 
guarantees a good correlation of Us and w^ and C^^ is of order 
one. 

For larger Rayleigh numbers the results approach the asymp- 
totic scaling 



(Rezon)p ~ (Repoi)p . 
When using Eq. (30) and Eq. (35) this suggests that 
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which is more or less confirmed by the correlation shown in 
Fig. 12. The small scale turbulent motions that develop for 
highly supercritical convection tend to destroy the correlation 
(Christensen, 2002; Showman et al., 2011). Once more, neither 
the ratio of (Rezon)p to (Repoi)p, nor the dependence of the cor- 
relation Cs(f) on <Ra*)p is significantly aff'ected by the diff'erent 
density stratifications employed. 

The loss of correlation explains why the ratio of total to 
non-zonal kinetic energy shown in Fig. 9c decreases for larger 
Rayleigh numbers. Our scaling of C^^ diff'ers from Showman 
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et al. (2011), who suggest an exponent of -1/2. This also ex- 
plains the difference in the scaHng of the Rossby number men- 
tioned above. This difference may come from specific assump- 
tions used in the model by (Kaspi et al., 2009). While these 
authors employ a more realistic equation of state, they also 
consider a simplified expression of the viscosity and neglect 
viscous heating. As Jones and Kuzanyan (2009) have demon- 
strated, viscous heating can contribute up to 50% to the heat 
budget, and is therefore likely to influence the average flow 
properties in the strongly nonlinear regime. 

4.2. Density stratification and vorticity production 

The z-component of the vorticity equation allows to analyze 
how the correlation between Us and w^ comes about. Taking the 
z-component of the curl of the momentum equation (5) leads to 
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where D/Dt denotes the substantial time derivative and co^ is 
the vorticity component along the axis of rotation. For small 
Ekman numbers, the first two terms on the right hand side dom- 
inate. In the Boussinesq approximation, the divergence of u 
vanishes and thus only the first term can contribute. This is 
called the vortex stretching term because it describes the vortic- 
ity changes in a fluid column experiencing the height gradient 
of the container (e.g. Schaeffer and Cardin, 2005; Glatzmaier 
et al., 2009). Because of the curved boundary, the absolute 
value of the gradient and thus the stretching effect increases 
with s, leading to the prograde tilt of the convective columns 
(e.g. Busse, 1983, 2002). At larger Rayleigh numbers, these 
columns loose their integrity along z and turbulent effects rather 
than the boundary curvature start to influence the flow. As re- 
ported above, this leads to the decrease in the correlation C^^. 

Besides this classical mechanism, compressibility brings a 
new vorticity source. As explained by Evonuk (2008) and 
Glatzmaier et al. (2009), rising (sinking) plumes generate nega- 
tive (positive) vorticity due to the compressional source (2E"^ -h 
6g)^)V • u in Eq. (39). In the 2-D anelastic simulations performed 
by Glatzmaier et al. (2009), the classical vortex stretching term 
is neglected and the new compressional term is the only source 
of vorticity. They report that the direction of the equatorial jet 
depends on the variations of the density background. According 
to Eq. (4), this new compressional source V-m = -(dlnp/dr) Ur 
is directly proportional to the density gradient. It is therefore 
expected to prevail closer to the outer boundary, where this gra- 
dient is strong. Moreover, since this mechanism is a local pro- 
cess, it may be less sensitive to the loss of correlation along 
a convective column with increasing Rayleigh number (Glatz- 
maier et al, 2009). 

At first sight, however, this is hard to reconcile with the fact 
that the zonal wind structure depends little on the density strat- 
ification in our simulations. Figures 11 and 12 have actually 
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Figure 13: Time-average of the correlation C^^ integrated over cylindrical sur- 
faces for simulations with Np = 10~^ (upper panel): Ra = 3 x 10^ (solid light- 
grey line), Ra = 5 X 10^ (dashed grey Hne), Ra = 8x10^ (dotted dark-grey line) 
and Ra = 1.5 x 10^ (dot-dashed black line); and for simulations with Np = 5 
(lower panel): Ra = 1 x 10^ (solid light-grey line), Ra = 1.8 x 10^ (dashed grey 
line), Ra = 4 X 10-^ (dotted dark- grey Hne) and Ra = 6 x 10^ (dot-dashed black 
line). The vertical line corresponds to the tangent cylinder. 



demonstrated that the decorrelation aff'ects all our cases to a 
similar degree, independently of the density background. 

In order to better understand the decorrelation. Fig. 13 com- 
pares the radial correlation profiles for nearly Boussinesq and 



Nn 



5 cases at various Rayleigh numbers. In the nearly 



Boussinesq cases, the correlation is predominantly lost in the 
interior but changes very little close to the outer surface, in 
agreement with the results of Christensen (2002). For the 
strongly stratified cases, the correlation indeed peaks close to 
the outer boundary as expected from the new compressional 
vorticity source, as long as the Rayleigh number remains mod- 
erate. However, on increasing the supercriticality, this particu- 
larly strong correlation is lost first and the maximum of C^^ is 
now located deeper in the interior. 

This is also illustrated by Fig. 14a, which shows a merid- 
ional cut of azimuthally averaged C^^. Figure 14b shows that 
the decorrelation goes along with a change in the vorticity pat- 
tern. In the inner part, the flow still shows a pronounced tilt that 
leads to Reynolds stresses and is clearly dominated by positive 
z-vorticity. Close to the outer boundary, convective motions 
are more radially oriented, without a preferred tilt or a clearly 
preferred sign. We attribute this to an increased mixing efl&- 
ciency in the very outer part of the flow (Aumou et al., 2007). 
The local turnover timescale of convection is a good proxy to 
quantify the mixing efficiency. We estimate this timescale to 
be Tto ~ Hp/ucom- This can be compared with a timescale as- 
sociated to zonal shear, defined by Tzf ~ ^/i^zonai (where S is 
the typical width of the jets). To yield a significant correlation 
between u, and u^, the shear timescale needs to be lower than 
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Figure 14: a) Correlation C^^ in the meridian plane, b) Vorticity along the axis of rotation cxj^ displayed in the equatorial plane. From a model with Np = 5, 
Ra = 4 X 10^. Positive values are rendered in red, negative ones in blue. The dashed lines in both panels correspond to the radius where the turnover timescale of 
convection becomes equal to the timescale associated with zonal shear (see Eq. 40). 

the turnover timescale, so that a fluid parcel could be diverted 
in azimuthal direction before loosing its identity, leading to 



Hn 



Tto > Tzf : 



^zonal 



(40) 



Figure 15 shows that this condition is fulfilled in the inner part 
of the layer but fails in a very thin layer close to the surface 
(approximately the outermost 5%). We have marked the ra- 
dius where both timescales become equal by a dashed circle in 
Fig. 14. 

As demonstrated above for stronger stratifications (see 
Fig. 5), Hp decreases outward while Wconv increases. This leads 
to a shorter turnover timescale close to the surface. Therefore, 
surface plumes that could potentially help to generate vorticity 
via the compressional source have a too short lifetime to sig- 
nificantly foster Reynolds stresses. In the deeper part, the extra 
eflficiency gained by compressible eff'ects is limited, explaining 
why zonal jets of compressible and Boussinesq simulations are 
found to be similar. 




2.2 2.3 

Radius 



5. Conclusion 

We have investigated the eff'ects of compressibility in 3-D 
rapidly rotating convection. Following Gilman and Glatzmaier 
(1981), Jones and Kuzanyan (2009) and Kaspi et al. (2009), 
all the simulations have been computed under the anelastic ap- 
proximation. As we focus on the eff'ects of the density stratifi- 
cation, we have deliberately chosen a moderate Ekman number 
of 10"^, which allowed us to reach more supercritical Rayleigh 



Figure 15: Radial profiles of the two timescales defined in Eq. (40): turnover 
timescale of convection (dashed black line) and timescale associated with shear 
(solid grey line). These profiles have been derived from an azimuthal average 
of the values in the equatorial plane for a model with Np = 5 and Ra = 4 x 10-^ 
(same model as in Fig. 14). 
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numbers. From the simulations computed in this parametric 
study, we highhght the following results: 

• Concerning the onset of convection, an increase of the 
density stratification is accompanied with a strong con- 
finement of the convective columns close to the outer 
boundary as well as an increase of the critical azimuthal 
wavenumber, in agreement with the results of the linear 
stability of Jones et al. (2009). This gradual outward shift 
of the onset of convection is explained by the radial in- 
crease of buoyancy. 

• When increasing supercriticality, the solutions which sim- 
ply drift close to onset, start to vacillate, then become 
intermittent and finally chaotic (Grote and Busse, 2001). 
While these diff'erent regimes are clearly separated in the 
Boussinesq cases, the transitions occupy a shorter fraction 
of the parameter space with increase of density stratifica- 
tion. 

• For stronger background density stratifications, the con- 
vective flow amplitude increases outward, while the 
lengthscale decreases. Concerning zonal winds, Boussi- 
nesq and anelastic simulations show very similar proper- 
ties, provided the convection is strongly driven: both the 
extent and the number of jets become fairly independent 
of the density background. This suggests that zonal flows 
follow a universal behaviour, independent of the density 
stratification, and may also explain why Boussinesq mod- 
els were already very successful in reproducing the mor- 
phology of the jets observed in gas giants (e.g. Heimpel 
et al, 2005). 

• In the strongly nonlinear regime, the solutions seem to 
follow the same asymptotic scaling laws, independently 
of the density stratification. Convective Rossby number 
and Nusselt number follow previously studied laws (Chris- 
tensen, 2002). The scaling of the total Rossby number 
diff'ers slightly from that suggested by Showman et al. 
(2011), which may be attributed to the diff'erent model 
used by these authors. The obtained scaling laws al- 
low us to extrapolate our simulations to Jupiter. Taking 
an internal heat flux of 5.44 W.m'^ (Hanel et al, 1981), 
Q = 1.75 X 10-4 s-^ Cp = I3x 10^ J.kg-^K"^ and d = 
1.2 X 10"^ kms, the evolution models of Guillot et al. (2004) 
suggest that the flux-based Rayleigh number decreases 
from Raj ~ 5 x 10"'^ at the surface to Raj ~ 3 x 10"^^ 
at the bottom of the molecular envelope (i.e. t] = 0.85). 
This leads to a mass-weighted flux based Rayleigh num- 
ber of (RaJ)p ~ 1.2 X 10"^^ Our scaling law (Eq. 30) then 
predicts a Rossby number of 6 x 10""^, which corresponds 
to an average zonal velocity around 1 m.s"^ In agreement 
with the previous extrapolations of Showman et al. (201 1), 
this value is much weaker (about one order of magnitude) 
than the surface zonal winds of gas giants. This rescaling 
is however tentative as the scaling laws may still depend 
on other parameters that have been kept constant in this 



study, for example Ekman number, Prandtl number or as- 
pect ratio, and further simulations would be required to 
clarify this. 

• The zonal jets are maintained by Reynolds stresses, which 
rely on the correlation between zonal (w^) and cylindri- 
cally radial (us) flow components. The gradual loss of 
this correlation with increasing supercriticality hampers all 
our simulations, independently of the background density 
stratification. However in the strongly stratified simula- 
tions, the correlation is first lost in the outer part of the 
flow, where compressibility efl'ects are most significant. 
Because of the short lengthscales and the large flow ampli- 
tudes found there, the fluid parcels loose their identity be- 
fore interacting efficiently with the zonal flows. This may 
explain why the additional compressional source of vortic- 
ity has not the important efl'ect expected by Evonuk (2008) 
and Glatzmaier et al. (2009) in our simulations. How- 
ever, the influence of the difl'erent vorticity contributions 
are hard to disentangle in 3-D compressible convection, 
stressing the need of future theoretical work to address this 
question. 

All the results obtained in this study have been derived in 
a non-conducting layer, where magnetic efl'ects have been ne- 
glected. Nevertheless, as suggested by Liu et al. (2008), the 
Ohmic dissipation associated with the outer semi-conducting 
region could strongly limit the extent of the surface jets. Ad- 
dressing this problem would require 3-D dynamo simulations 
of compressible convection with variable conductivity (Stanley 
and Glatzmaier, 2010; Heimpel and Gomez Perez, 2011). 
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Table 2: List of numerical simulations presented in this work. For all the simulations, E = 10 and Pr = 1. 



Ra/Racrit 


<Ro>p 


<Ropoi>p 


Kf/E^^n 


Nu 




N 


p = 10 


r2, symbol=red circle 




1.0006 


4.49 X 


10-5 


2.49 X 


10-5 


0.002 


1.0002 


1.0063 


1.03 X 


10-4 


5.71 X 


10-5 


0.011 


1.0009 


1.035 


2.35 X 


10-4 


1.26 X 


10-4 


0.058 


1.0045 


1.15 


5.04 X 


10-4 


2.45 X 


10-4 


0.231 


1.016 


1.72 


1.25 X 


10-3 


4.41 X 


10-4 


0.586 


1.045 


2.01 


1.52 X 


10-3 


4.96 X 


10-4 


0.640 


1.053 


2.30 


1.98 X 


10-3 


5.86 X 


10-4 


0.700 


1.067 


2.87 


2.85 X 


10-3 


7.19 X 


10-4 


0.773 


1.087 


3.45 


3.37 X 


10-3 


8.23 X 


10-4 


0.788 


1.104 


4.60 


5.05 X 


10-3 


1.24 X 


10-3 


0.816 


1.1800 


5.75 


7.18 X 


10-3 


1.58 X 


10-3 


0.843 


1.23 


7.47 


1.16X 


10-2 


2.03 X 


10-3 


0.902 


1.31 


8.62 


1.44 X 


10-2 


2.35 X 


10-3 


0.916 


1.37 


10.35 


1.95 X 


10-2 


2.97 X 


10-3 


0.929 


1.51 


11.50 


2.24 X 


10-2 


3.38 X 


10-3 


0.933 


1.61 


14.37 


3.09 X 


10-2 


4.51 X 


10-3 


0.942 


1.89 


17.25 


3.90 X 


10-2 


5.62 X 


10-3 


0.945 


2.20 


20.12 


4.51 X 


10-2 


6.67 X 


10-3 


0.943 


2.52 


23.00 


5.20 X 


10-2 


7.96 X 


10-3 


0.939 


2.89 


28.75 


6.34 X 


10-2 


1.04 X 


10-2 


0.932 


3.65 


34.50 


7.30 X 


10-2 


1.27 X 


10-2 


0.924 


4.44 


40.25 


8.30 X 


10-2 


1.52 X 


10-2 


0.915 


5.29 


46.00 


9.10 X 


10-2 


1.77 X 


10-2 


0.903 


6.04 


51.75 


9.98 X 


10-2 


2.00 X 


10-2 


0.893 


6.81 


57.50 


1.09 X 


10-1 


2.24 X 


10-2 


0.883 


7.60 


69.00 


1.21 X 


10-1 


2.73 X 


10-2 
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9.22 


86.25 


1.42 X 


10-1 


3.40 X 


10-2 
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11.53 
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1.79 X 


10-1 


4.78 X 


10-2 


0.845 
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2.03 X 


10-1 


5.74 X 


10-2 
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19.50 




N 


^ — 7 


symbol=l 


blue triangle 
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10-5 


2.63 X 


10-5 


0.001 


1.0001 


1.0035 
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10-5 
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10-5 
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10-4 
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3.94 X 
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<Ropoi>p 



V^kin 



Nu 



A/p = 1, symbol=green triangle 



1.0010 
1.0068 
1.0145 
1.024 
1.062 
1.15 
1.35 
1.44 
1.54 
1.73 
1.93 
2.51 
2.89 
3.28 
3.86 
4.83 
5.79 
7.72 
9.66 
11.59 
13.52 
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19.32 
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28.98 
38.64 
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1.0020 
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1.0072 
1.013 
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1.63 
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32.70 
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1.0038 
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1.95 X 10-4 
2.54 X 10-4 
4.23 X 10-4 
6.99 X 10-4 
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1.48 X 10-3 
1.78 X 10-3 
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3.11x10-3 
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6.74 X 10-2 
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2.71 X 10-3 
3.57 X 10-3 
4.43 X 10-3 
6.28 X 10-3 
8.39 X 10-3 
1.08 X 10-2 
1.31 X 10-2 
1.53 X 10-2 
2.02 X 10-2 
2.48 X 10-2 
3.04 X 10-2 
3.95 X 10-2 
4.87 X 10-2 



0.002 
0.014 
0.033 
0.055 
0.156 
0.357 
0.537 
0.619 
0.674 
0.748 
0.795 
0.858 
0.885 
0.905 
0.925 
0.922 
0.931 
0.937 
0.927 
0.916 
0.909 
0.904 
0.896 
0.893 
0.875 
0.873 
0.871 



Np - 3, symbol=magenta diamond 



4.24 X 10-5 
6.35 X 10-5 
9.32 X 10-5 
1.27 X 10-4 

2.29 X 10-4 
3.14x10-4 
5.17x10-4 
1.75 X 10-3 
2.91 X 10-3 
3.80 X 10-3 

6.30 X 10-3 
9.09 X 10-3 

1.47 X 10-2 
2.07 X 10-2 
2.62 X 10-2 
3.21 X 10-2 
3.77 X 10-2 
4.91 X 10-2 
6.66 X 10-2 

8.48 X 10-2 
1.03 X 10-1 

1.25 X 10-1 
1.47 X 10-1 
1.70x10-1 



3.32 X 10-5 
5.15x10-5 
7.21 X 10-5 
1.01 X 10-4 
1.75 X 10-4 
2.30 X 10-4 
3.35 X 10-4 
7.85 X 10-4 
1.09x10-3 

1.27 X 10-3 
1.75 X 10-3 
2.15x10-3 
2.99 X 10-3 
3.95 X 10-3 

4.90 X 10-3 
5.84 X 10-3 

6.91 X 10-3 
8.94 X 10-3 

1.28 X 10-2 
1.74x10-2 
2.17x10-2 
2.70 X 10-2 
3.28 X 10-2 
4.09 X 10-2 



0.001 
0.002 
0.004 
0.008 
0.037 
0.111 
0.307 
0.675 
0.770 
0.812 
0.866 
0.899 
0.920 
0.924 
0.922 
0.924 
0.922 
0.922 
0.916 
0.908 
0.905 
0.902 
0.900 
0.886 



Np - 5, symbol = orange square 



2.49 X 10-5 
3.56 X 10-5 
9.09 X 10-5 
1.50x10-4 
2.71 X 10-4 
4.61 X 10-4 
8.20 X 10-4 
1.32x10-3 
1.96x10-3 
2.82 X 10-3 
3.92 X 10-3 
5.04 X 10-3 
9.65 X 10-3 
1.30x10-2 
1.82 X 10-2 
3.09 X 10-2 
3.99 X 10-2 
5.03 X 10-2 
6.26 X 10-2 
6.94 X 10-2 
7.85 X 10-2 



2.24 X 10-5 
3.21 X 10-5 
8.06 X 10-5 
1.33 X 10-4 
2.19x10-4 
3.10x10-4 
4.50 X 10-4 
6.19x10-4 
8.10x10-4 
1.08 X 10-3 
1.44x10-3 
1.77 X 10-3 
2.92 X 10-3 
3.56 X 10-3 
4.45 X 10-3 
6.65 X 10-3 
8.17x10-3 
1.01 X 10-2 
1.26x10-2 
1.40x10-2 
1.63 X 10-2 



0.000 
0.001 
0.003 
0.013 
0.124 
0.351 
0.548 
0.661 
0.726 
0.751 
0.759 
0.766 
0.821 
0.853 
0.881 
0.904 
0.911 
0.914 
0.914 
0.915 
0.911 



1.0001 
1.0008 
1.0018 
1.0029 
1.007 
1.009 
1.028 
1.036 
1.044 
1.058 
1.07 
1.14 
1.23 
1.26 
1.37 
1.49 
1.69 
2.17 
2.79 
3.42 
4.14 
4.86 
6.41 
7.84 
9.69 
12.86 
15.97 



1.0003 

1.0007 

1.0013 

1.0025 

1.0076 

1.0134 

1.024 

1.074 

1.11 

1.13 

1.20 

1.27 

1.45 

1.76 

2.09 

2.42 

2.81 

3.51 

4.80 

6.30 

7.84 

9.47 

11.76 

14.14 



1.0005 
1.0011 
1.0066 
1.017 
1.041 
1.061 
1.104 
1.163 
1.24 
1.37 
1.57 
1.77 
2.37 
2.67 
3.06 
4.01 
4.68 
5.45 
6.52 
6.87 
8.04 
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